Molecular orbital calculations of two-electron states for P donor solid-state spin qubits 
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We theoretically study the Hilbert space structure of two neighbouring P donor electrons in 
silicon-based quantum computer architectures. To use electron spins as qubits, a crucial condition 
is the isolation of the electron spins from their environment, including the electronic orbital degrees 
of freedom. We provide detailed electronic structure calculations of both the single donor electron 
wave function and the two-electron pair wave function. We adopted a molecular orbital method for 
the two-electron problem, forming a basis with the calculated single donor electron orbitals. Our 
two-electron basis contains many singlet and triplet orbital excited states, in addition to the two 
simple ground state singlet and triplet orbitals usually used in the Heitler-London approximation to 
describe the two-electron donor pair wave function. We determined the excitation spectrum of the 
two-donor system, and study its dependence on strain, lattice position and inter donor separation. 
This allows us to determine how isolated the ground state singlet and triplet orbitals are from the 
rest of the excited state Hilbert space. In addition to calculating the energy spectrum, we are also 
able to evaluate the exchange coupling between the two donor electrons, and the double occupancy 
probability that both electrons will reside on the same P donor. These two quantities are very 
important for logical operations in solid-state quantum computing devices, as a large exchange 
coupling achieves faster gating times, whilst the magnitude of the double occupancy probability can 
affect the error rate. 

PACS numbers: 03.67.Lx, 71.55.Cn, 85.30.De 



I. INTRODUCTION 

Recently several designs for silicon-based quantum 
computer architectures have been proposediSi 3 *^*^ In 
this work we concentrate our efforts on the Kane model, 1 
which exploits a qubit array of nuclear spins of 31 P 
dopants embedded within a silicon crystal matrix. The 
model is based on the use of 31 P nuclear spins as qubits, 
with the donor electrons functioning to mediate control of 
single qubit operations via the hyperfine interaction, and 
interaction between individual qubits via the exchange 
interaction, and permit read-out of nuclear spin states. 

Performing logical operations on either electron-spin 
or nuclear-spin solid-state qubits requires precise con- 
trol over single and two-qubit unitary operations, which 
corresponds to precise control over the electron-electron 
exchange interaction and the electron-nucleus hyperfine 
interaction in the Kane quantum computer. Here we cal- 
culate the exchange interaction as a function of the two 
donors' relative positions in the lattice and strain. We use 
multi-valley effective mass theory to calculate the single 
donor electron wave functions, these single donor orbitals 
combine to form our two-electron basis. This theory in- 
corporates the Si crystal lattice effects by including the 
Si crystal Bloch functions into our single donor electron 
basis. Instead of using the Heitler-London (H-L) approx- 
imation, which has been used extensively in the literature 
so far for impurities in SijSiiLiSiiAiiSi we describe our two- 
donor system using a rigorous molecular orbital method 
which employs our multi-valley single donor orbitals to 
form our two-electron basis, to calculate the exchange 
coupling more accurately. 



An important feature necessary for quantum comput- 
ing is to have well-characterized qubits, and for the two 
qubit case this is the ground state singlet and triplet two- 
electron states. It is meaningful to study the degree of 
proximity these targeted ground state orbitals are to the 
rest of the unwanted excited state Hilbert spaced This 
energy separation gives us an estimate for the conditions 
under which adiabaticity can be attained. We have pur- 
sued this goal using a molecular orbital method, which 
enables us to calculate a large number of two-electron 
energy levels (144), in the energy spectrum for our two 
donor system. 

We used a molecular orbital method which includes 
the single donor ground state and first five excited states 
at each donor to form the two-electron basis. This yields 
a basis of 78 singlet states and 66 triplet states. This 
method not only gives us the exchange coupling, which 
is the difference between the ground singlet and triplet 
two-electron states, but also the spectrum of energy lev- 
els for the two donor system. For comparison we also 
calculated the H-L exchange coupling using just the sym- 
metrized and anti-symmetrized products of the single 
donor ground states, and the Hund-Mulliken (H-M) ex- 
change coupling which in addition to the H-L states, 
includes the two doubly occupied single donor ground 
states at each donor, in our two-electron basis. We cal- 
culated the exchange coupling and energy spectrum for 
our two-electron system as a function of donor position 
and strain. In addition, we also calculated the prob- 
ability that both electrons will be on the same donor. 
This is also an important parameter for quantum gate 
operations, as these doubly occupied states can become 
a potential source of error. Several authors have similarly 
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studied the exchange coupling, double-occupancy errors 
and adiabaticity of spin qubits in a number of different 
solid-state quantum computing architectures 

Much attention has been devoted to modeling 
the hyperfine and exchange interactions in these 
devicesi 9iP' 1 ?i 11 ' 12 i 1 7i 1 9i 1 ?i 2l ' ) i 21 Inter valley interference be- 
tween degenerate conduction band minima in Si has 
been shown to lead to oscillations in the exchange cou- 
pling as a function of the donor pair positioning in the 
latticeAifiiiiiiS This poses serious problems for the fab- 
rication of these devices, and leads to an extreme sen- 
sitivity of the exchange energy on the relative orienta- 
tion of the P atoms. Koiller et alM* demonstrated that 
the introduction of external strain on the Si lattice par- 
tially lifts the valley degeneracy in the bulk Si. They 
showed that the inter valley effects could be reduced in 
some cases depending on the relative orientations of the 
donor pairs, whilst in other cases the donor exchange cou- 
pling remains oscillatory. The molecular orbital method 
we employ not only improves the calculation of the two 
donor electron wave functions, but we also use a more 
flexible basis than previous studie a?' 1 ^ 11 ' 12 to calculate 
the single donor electron wave functions, which are used 
to construct the two-electron basis. 



II. QUANTUM CHEMICAL MODELS 

We advance beyond the simple H-L model for the 
two-electron wave function in the Kane device that has 
been previously consideredA^iiSiiiiiS In the molecular 
orbital method we use the single donor wave functions 
to form our basis states, and solve the 6-D Schrodinger 
equation for the two electrons through a direct matrix 
diagonalisationii^ 

In the simplest case, the H-L approximation, the 
donor pair wave function is modeled as the symmetrized 
and anti-symmetrized products of the two single donor 
ground state wave functions ("Ai" states) at each P nu- 
cleus, to form our singlet and triplet states respectively. 
In the H-M approximation, in addition to the two H- 
L states, the H-M basis incorporates the two "ionised" 
or "polarised" doubly occupied ground states, at each 
donor. 



For the molecular orbital calculation we extended these 
bases to also include the first five excited states for each 
donor in our basis, in addition to the single donor ground 
states. This was chosen so that our basis included the six 
symmetry ground states for the single P donor, {A\^2 
and E states). We performed calculations for the ex- 
change coupling to see the effect that this larger basis 
has on lowering the energy of both the singlet and triplet 
ground states, and to improve upon and test the validity 
of using H-L theory to model the two-donor system over 
a range of device parameters. 

We get a basis for our two-electron system which con- 
sists of the spatially symmetric singlet states, and anti- 
symmetric triplet states. Because the spin part of the 
singlet and triplet states are orthogonal, we can consider 
the singlet and triplet bases independently. Using six 
single donor orbitals on both qubits, we can form 78 
singlet states and 66 triplet states. The molecular or- 
bital method has advantages over some quantum chem- 
ical methods as it includes the correlation between the 
two electrons, by virtue of including many two-electron 
orbitals to minimise the energy of the system. 
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FIG. 1: Co-ordinate geometry of our two-electron problem. 

We show the geometry of our two-electron problem in 
Fig. ^for two P donors, Q\ and Q2, embedded in the Si 
lattice, (the origin is at Qi). For the singlet (symmetric 
spatial orbitals) basis we form the following two-electron 
wave functions, from our basis of single donor orbitals, 
^ e Q i {r) and ^ e Q 2 {r — R), (e n th state at Qi, and e m th 
state at Q2 respectively): 



*1-21 



*f 3 -78 



for n — to 5 and m = n to 5. 

V2(1 = Snm) - *)*&(^ -r)+ *q->2 - mn m M - m 

for n — to 5 and m = n to 5. 

1 (ri)*Sr(ra — R) + (r 2 )*o m (ri - R) 

v/2(l + |S nm | 2 ) L QlV QaV Ql1 * l 

for n — to 5 and m — to 5, 



3 



where S nm = J dr^r^r - R). 



Here we see that the two-electron singlet donor wave 
functions 5 , i_ 2 i and 'I' 22 -42; are the doubly occupied sin- 
glet states located at Q\ and Q2 respectively. The two- 
electron states, V I , 43_ 78 , are the "Heitler-London like" sin- 
glet states formed from the single donor ground state and 



excited state wave functioins. 

Similarly for the triplet (anti-symmetric spatial or- 
bitals) basis we obtain the following two-electron wave 
functions: 



*1-16 
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*31-66 



1 

for 
1 



*Q"(r-i)^(^)-^(^)*S:(^). . 

n = to 5 and to — n + 1 to 5, (to 7^ n). 

1 to 5, (to 7^ n). 



to = n + 



72 

for n = to 5 and 

V2(l-|5„ m | 2 ) L W 1J * V Ql v 2; Q2 v 1 

for n = to 5 and to = to 5. 
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It is clear that the singlet and triplet bases contain 
the original H-L states, *f?f 3 in the singlet basis, and 'f^ 
in the triplet basis. For the H-M calculation we include 
the two additional "ionised" or doubly occupied ground 
states, Vtf and ^f 2 in our singlet basis. In the extended 
molecular orbital basis, we consider all 78 singlet states 
and 66 triplet states in our two donor electron Hamilto- 
nians for the singlet and triplet bases respectively. 

III. NUMERICAL METHOD 
A. Solution of the single donor wave function 

To obtain our two-electron states we first need to eval- 
uate the single donor wave functions at each donor to 
use in the two-electron basis. We did this in the case of 
no strain and with uniaxially strained Si. We calculated 
the single donor orbitals using multi-valley effective mass 
theory. We use a basis for the multi-valley single donor 
wave function which includes the full Bloch structure at 
each conduction band minimum, in our basis functions. 

The Kohn-Luttinger— form of the wave function for a 
donor (where we define the P nucleus to be at the origin 
on a substituted Si atom site) is given by: 

*(r) = ^^(^(r), 

6 

= ^^(rK^^ + Ro), (1) 



where we choose Ro = (a°/8)(l, 1, 1) = i?g, 23 ' 2A as we 
take the origin to be at an substitutional P donor site, 
here a = 0.543nm is the length of the unit cell. The term 
Rq arises because Si has two atoms per lattice point in 
the unit cell, where the Si atoms are chosen to be dis- 
placed a distance of ±i?Q away from each lattice point in 
the unit celli^iSi For substitutional donors Ro = ±i?g. 

F^(r) is the donor envelope function at the conduc- 
tion band minimum, fc^, and ip® (r) is the silicon crystal 
Bloch function at the conduction band 
where fci = ±k z = (0, 0, ±k)2ir/a° etc. and k = 0.85. 

The multi-valley effective mass equation for a P donor 
in Si under strain is>2£ 

6 

a^-^-TO-iV) + U(r) + H stram - E] 

xfW(r)=0, (2) 

where: 

cT\ _cP_ 

1 \dx 2 dy 2 ) ^ dz 2 - 1 

= T 2 (-tV), 

\dx 2 dz 2 ) ^ dy 21 
= T 4 (-iV), etc. 

Here T M are the anisotropic kinetic energy terms, due to 
the anisotropy of the conduction band minima in Si. The 
impurity potential, U(r), is the potential term due to the 
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effective +1 charge of the P nucleus in the Si lattice. Here 
we model the impurity potential as a screened Coulombic 
potential, U(r) =2/r. We are using atomic units, where 
the unit of length [as] = h 2 e/m±e 2 = 31.667A and unit 
of energy [E B ] = m_ L e' i /2h 2 e 2 = 19.9436meV, where 
e = 11.4a.u. and 7 = mx/m\\ = 0.2079ii£ H stra in is 
the potential due to uniaxial strain along the z-direction 
which we will define later. 

We expanded the donor electron envelope wave func- 
tion, fM, in a basis of the single- valley zero field en- 
velope wave functions, Fj , at each of the six conduc- 
tion band minimum, k M . Here Fj are the eigenfunc- 

tions of the single-valley zero field Hamiltonian, = 
Tfj,{— iV) + U(r). We have discussed previously^ how we 

obtained the single- valley zero-field wave functions, Fj , 
by expanding these single- valley wave functions in a basis 
of deformed hydrogenic orbitals. 

In Eq. we expanded the donor electron wave func- 
tion, ^(r), in a basis of the donor electron envelope func- 
tions, F^\ at each minimum. In addition we can also 
expand F^' in our basis of single-valley donor electron 
wave functions, Fj : 

6 

= E e!k, '- r ^(r + Ro)E C, i Al) ^ ) «' 

(3) 

where ^(r) = £ B^ m ^ m (x, y, z, a,(3), 

n,l,m 

and C- are the expansion co-efficients for our ba- 
sis functions Fj(r). We see that the single- valley 

envelope functions, Fj, are in turn, also a sum 
of basis functions: the deformed hydrogenic orbitals, 

$nlm( x ' y> z ' a "^)' gi ven already in a previous paper>i£ 
The co-efficients B^^ are determined already since 
Fj^ (r) are the eigenf unctions of the single- valley Hamil- 
tonians, ie. H^F^(r) = E°F^\r). 

Note that including the expansion co- efficients m 
Eq. © is a generalisation of the calculations of Wellard 
et al&iSL where we have removed the restriction that the 
donor wave function be composed of equal contributions 
from the six conduction band minima. Clearly this re- 
striction breaks down when an external strain is applied 
as this will break the degeneracy of the six conduction 
band minima. 

So now Eq. (J2J becomes: 

£ e^-M- ]T Cf [HM + H stram - E]F^ (r) = 0. 

(4) 



We now multiply Eq. (J2J by F* (r) and integrate over 
r. The orthonormality of this basis is enforced by the 
e i(k M -k, y ).r terms w hi cn appear in the matrix elements, 
and due to their rapidly oscillating nature average to zero 
unless k M = 

In the standard effective mass treatment, the inter val- 
ley mixing terms which couple the envelope functions at 
different conduction band minima in the above approxi- 
mation are neglected, and six independent equations are 
obtained. For the higher donor excited states this is a 
valid approximation, as their energies agree quite well 
with calculations using only single valley effective mass 
theory*^ However, we need to consider the inter valley 
coupling for the donor ground state, which has the ef- 
fect of lifting the six- fold degeneracy of the IS states 
predicted by the one-valley effective mass equations. In 
order to obtain the correct symmetry states for the donor 
ground state, of a singlet (A\), a triplet (T2) and doublet 
(E), we add empirically determined parameters to our 
Hamiltonian, as was done by Koiller et a/^i 

Hence the multi valley effective mass Eq. (0} becomes: 

+ / dvF; {u \v)[H stram ]F^\v) 

j J 

= E5^5 13 C[ V) . (5) 

{0, if v = v., 

— 2.1934meV, if [i, v are on perpendicular 
symmetry valleys, 
— 1.535meV, if fi, v are on opposite 
symmetry valleys. 
Here we also scale the single valley ground state energy, 
E\ = — 35.19meV, to reproduce accurately the experi- 
mental splitting for the P donor ground state. We only 
scaled the single valley ground state energy because the 
single valley calculation reproduces the higher excited 
state energies reasonably accurately^ 

When we consider the effect of strain on the single 
donor orbitals, we consider its effect only on the lowest 
six energy states, because later we construct our two- 
electron basis from these lowest six single donor states. 
Here we follow the treatment of Koiller et alm^ and in- 
troduce the relative energy shifts due to uniaxial strain 
along the z-direction, in terms of a dimensionless val- 
ley strain parameter, \- I n their paper they discuss the 
physical relevance and tuning of this parameter. For our 
purposes we consider four cases of the strain parameter, 
corresponding to % = 0, —1, —5 and —20. Negative val- 
ues of x correspond to tensile strain, which favours the 
z-envelopes energetically, and % = — 20 represents the 
realistic situation of Si grown over relaxed Sio.sGeo^^ 
To evaluate the H stra in terms we first need to define 
ji = 1, 2, 3, 4, 5, 6 to correspond to the z, —z, y, —y, x, —x 
valleys respectively. Now the H str ain terms in Eq. (J5J) 
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FIG. 2: Contour plot of the ground state electron density in 
the yz-plane for A\ state without any strain applied. Here 
the P nucleus is located at the origin. 




FIG. 3: Contour plot of the ground state electron density in 
the yz-plane with a strain parameter \ — —20. Here the P 
nucleus is located at the origin. 

become: 

6 . 

S V E C i / [H stram ] ff>(r) 

where Ac = 2.16meV is used to be consistent with 
Koiller et alJ± 

We plot the ground state electron density without any 
external strain, x = 0, and with strain applied, \ — ~20, 
in Fig. to ^ using Eq. J2J. Here the Bloch func- 
tions are obtained using the empirical pseudo-potential 
technique^SI and the P donor envelope functions are 
obtained from the multi- valley effective mass equations. 
Figure [21 is a contour plot of the ground state electron 
density in the yz-plane for the symmetric A\ state corre- 
sponding to zero strain, where the contribution from all 
six valleys are equivalent. However with a strain applied 
in the z-direction, we see in Fig. |3 where we plotted the 
electron density in the yz-plane, the effective Bohr ra- 



FIG. 4: Contour plot of the ground state electron density in 
the zy-plane with a strain parameter \ — —20. Here the P 
nucleus is located at the origin. 

dius in the z-direction is reduced. This is because with 
an external strain applied, the six-valley degeneracy of 
the symmetric A\ ground state is broken. This can be 
seen from Eq. (jfJJl, and the lowest energy state is the one 
in which the effective Bohr radius in the direction par- 
allel to the strain is reduced, ie. the F^ Z (1S) states. 
In contrast in Fig. where we plotted the ground state 
density in the xy-plane, we see that the strain (applied 
in the z-direction) is equivalent in these two directions, 
and the effective Bohr radii has increased. 

Table[I]reports the energy splitting between the ground 
state and first excited state for a single donor electron 
for different magnitudes of strain applied. The energy 
levels become closer together when a strain is applied, 
and the ground u A\ n state is no longer degenerate in the 
six valleys, and we find that the F± z valleys become more 
favored. This leads to a smaller effective Bohr radius in 
the z-direction, and larger effective Bohr radii in the x, y- 
directions, which was demonstrated already in Fig.|Sjand 

El 

TABLE I: Energy splitting between the ground state and the 
first excited state for a single donor electron. 



X 


AE (meV) 





11.847 


-1 


8.316 


-5 


4.383 


-20 


3.378 



B. Solution of the two-electron donor pair wave 
function 

Once the single donor orbitals are known, we are then 
able to evaluate the 6-D two electron Hamiltonian matri- 
ces for both our singlet and triplet bases, H^e an d ^2e) 
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and the singlet and triplet overlap matrices, S s and S T . Here the 6-D two electron Hamiltonian operator is: 



H2e = 



,(r 2 )- 



2 



|ri-R| |r 2 -R| 



ri - r 2 



strain 



r 



Thus we need to evaluate both the singlet and triplet 
Hamiltonian matrix elements, { \I/ ^_ 78 1 | ''I' f _ 78 ) and 
H^el^T-m) respectively, and singlet and triplet 



overlap matrix elements, S(i,j) 
inter donor separation and strain. 



(^il^j), for varying 



Since our basis functions ^'q 1 (r) and (r — R) are 
eigenf unctions of the single electron Hamiltonian opera- 
tor at each donor, we use this property in evaluating the 
two electron Hamiltonian matrix elements: 



-V ams (r) - - + H stram (r) 



r-R 



H~~ H strain\*) 



*&(r-J*) 
I 



E%*%(r-R). 



Here when we calculate the matrix elements for the 
singlet and triplet Hamiltonian and overlap matrices, we 
retain the single plane wave part, e M' r , of the Bloch 
functions at each minima in the expansion for the single 
donor electron wave functions, in the integrands involv- 
ing these wave functions. This leads to the inherent os- 
cillations in the exchange energy due to the inter valley 
interference between these terms at the degenerate con- 
duction band minima. We still neglect the periodic part 
of the Bloch functions, itfc (r + Ro), in the integrands. 
It has been shown^ that this is an excellent approxima- 
tion, and it was impossible to distinguish between the 
results for the exchange coupling using this approxima- 
tion, and those including the detailed Bloch structure. 
We did this to make the calculations more tractable over 
a larger range of device parameters. 

Once we derived the matrix elements for both the sin- 
glet and triplet Hamiltonian and overlap matrices, we 
needed to solve a generalised eigenvalue problem for both 
the singlet and triplet case. This is because the two- 
electron states are not necessarily orthogonal, since the 
single electron wave functions, (r) and VPg™ (r — R), 
are not orthogonal. We have: 



H 2e c = ESc. 



(7) 



Here c is a vector of the coefficients of the two-electron 
basis functions. To solve this we first need to compute 
the Cholesky factorisation for the overlap matrix S, to 
give S — LL + . We did this using a standard numerical 
subroutine. Once we had obtained the Cholesky factori- 
sation, we used this to transform Eq. J7J into the stan- 



dard eigenvalue problem using another subroutine: 

[L~ 1 H 2e (L + )^ 1 ][L + c] = E[L+c\. (8) 

Once we derived the standard eigenvalue problem, we 
used a standard eigenvalue solver, to diagonalise Eq. (JSJ) 
to obtain the energies E, for the singlet and triplet states. 

The most computationally expensive task in our molec- 
ular orbital calculations is the computation of the 6-D 
two-electron integrals in the singlet and triplet Hamilto- 
nian matrix elements. In the singlet basis this required 
3081 6-D integrals to be performed, and for the triplet 
basis 2211 6-D integrals to be performed. We have re- 
duced this task greatly by only calculating the identical 
6-D integrals in both the singlet and triplet bases once. 
This means calculating 4131 6-D integrals in total. This 
is also a better numerical practice as it means that when 
we calculate the energy splitting between the ground sin- 
glet and triplet states, (which can be very small), we 
are using the same integral evaluations to calculate both 
quantities, Et and Es- Thus the exchange energy calcu- 
lated J = Et — Eg will be more accurate, as the same 
numerical errors will be involved in both quantities. 

We have also increased the efficiency and speed of 
our code by modifying a standard Monte Carlo subrou- 
tine used to numerically evaluate the 3-D and 6-D inte- 
grals. We did this because the integrals all require eval- 
uations of our single donor basis functions, ^g" (r) and 
^q™(r — R). We have greatly reduced the complexity 
and computing time for these calculations by evaluating 
these common basis functions on a grid, before inputting 
these functions into the Monte Carlo subroutine. This 
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FIG. 5: We compare the exchange coupling at lattice sites 
along the [010] or y direction for small magnitudes of R, cal- 
culated using our three quantum chemical models in (a): us- 
ing the H-L states, H-M states and our extended molecular 
orbital basis, for zero strain. In (b) we plot the singlet and 
triplet two-electron energy levels using our extended basis. 
Here we only consider values of R such that both P donors 
are on substitutional donor sites. 
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has provided a speed-up of our calculations of the order 
of 100 times. 
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IV. RESULTS AND DISCUSSION 
A. Results using full molecular orbital calculation 

We give the results for our three quantum chemical 
models for the two-electron states: the two H-L states, 
the four H-M states and our extended molecular orbital 
basis. Figures 02a) to|5Ja) show a comparison of the ex- 
change coupling obtained using our three methods, for 
varying inter donor separations, and for Q\ and Q2 lo- 
cated at lattice sites. The reason why we study the ex- 
change coupling in more detail for R greater than 14nm in 
Fig. to[|Ja) is because separations of about at least 
14nm are envisioned to be needed in order for metal- 
lic gates to be placed on top of and between adjacent 
qubits (currently the smallest width of the metallic gates 
that can be fabricated is about 10nm)i These gates pro- 
vide additional tuning of the electron density and P nu- 



FIG. 6: We compare the exchange coupling at lattice sites 
along the [010] or y direction, calculated using our three quan- 
tum chemical models in (a) for R > 14nm: using the H-L 
states, H-M states and our extended molecular orbital basis, 
for zero strain. In (b) we plot the ground state singlet and 
triplet energy separately for clarity, and in (c) we plot the rest 
of the excited two-electron energy levels using our extended 
basis. Here we only consider values of R such that both P 
donors are on substitutional donor sites. 

clear spin via the application of varying voltages to them. 
Here we consider the inter donor separations along the 
y or [010] direction only (see Fig. [TJ. This is because 
executing the full molecular orbital calculation is very 
computationally expensive. In the next section we use 
the H-M method to calculate the exchange coupling for 
many different orientations of Q\ and Q2 in the lattice. 

We can observe that as R increases the H-L calculation 
is more accurate as the two donors become further sep- 



8 



> 

a 



0.03 r 

0.025 

0.02 - 
0.015 

0.01 
0.005 - 
- 



» (a) 



14 



14.5 



15.5 16 
R(nm) 



16.5 



17.5 



All states in basis O 
Hund-Mulliken states X 
Heitler-London states + 



-99.8 
-100 
-100.2 
-100.4 
-100.6 
-100.8 
-101 
-101.2 
-101.4 
-101.6 



- (b) 



14.5 



15.5 16 
R (nm) 



16.5 



17.5 



Ground state singlet energy 
Ground state triplet energy 



(c) 



_90 -93.22 
-93.24 

-92 - 



15 16 17 

R(nm) 

Excited state singlet energy levels 
Excited state triplet energy levels 



> 
1 



0.04 p| 
0.035 - 

0.03 - 
0.025 

0.02 - 
0.015 

0.01 

0.005 

— 
14 



(a) 



14.5 



15.5 16 
R(nm) 



16.5 



17.5 



All states in basis O 
Hund-Mulliken states X 
Heitler-London states + 



-127 



-127.5 



-128.5 



(b) 



14.5 15 



15.5 16 
R (nm) 



16.5 



17.5 



Ground state singlet energy 
Ground state triplet energy 





-90 r 




-95 




-100 


> 


-105 


£ 






-110 


bo 






-115 


a 






-120 




-125 




-130 



... (c) 



-123.8 



-123.9 



15 



16 

R(nm) 

Singlet energy levels 
Triplet energy levels 



FIG. 7: Comparison of the exchange coupling for \ — — 1 in 
(a) along the [010] or y direction. In (b) we plot the ground 
state and triplet state energies separately for clarity, and in 
(c) we plot the rest of the excited two-electron energy levels 
using our extended basis. Here we only consider values of it 
such that both P donors are on substitutional donor sites. 



FIG. 8: Comparison of the exchange coupling for \ = —5 in 
(a) along the [010] or y direction. In (b) we plot the ground 
state and triplet state energies separately for clarity, and in 
(c) we plot all the singlet and triplet two-electron energy levels 
using our extended basis. Here we only consider values of R 
such that both P donors are on substitutional donor sites. 



arated, and it becomes a better approximation to treat 
the two donors as a superposition of the single electron 
ground state wave functions centered at each donor. We 
demonstrate this in Fig.|SJa) and^a) using smaller inter 
donor separations (R < 12nm), and larger inter donor 
separations (R > 12nm) respectively, for \ — 0. We 
found that the exchange coupling is improved substan- 
tially for the smaller inter donor separations using our 
full molecular orbital calculations. We can also see that 
even if we just include the "doubly occupied" states in 
our H-M calculation we get a significant improvement in 



the exchange energy over H-L theory, when we compare 
it with the full molecular orbital calculation. 

In part (b) of Fig. [S]to[!2| we show the exchange split- 
ting between the ground singlet and triplet states using 
the full molecular orbital calculations. The results for 
the two P donor-pair wavefunction give a two-electron 
ground state of the order -98 to -lOOmcV for R > 14nm 
and x = 0, (shown in Fig. HJb)). The energy of two iso- 
lated P atoms should be in the order of -91meV. When 
we calculate the singlet and triplet energy levels, the 6-D 
and 3-D integrals involve both repulsive, direct Coulomb 
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FIG. 9: Comparison of the exchange coupling for x = —20 in 
(a) along the [010] or y direction. In (b) we plot the ground 
state and triplet state energies separately for clarity, and in 
(c) we plot all the singlet and triplet two-electron energy levels 
using our extended basis. Here we only consider values of R 
such that both P donors are on substitutional donor sites. 



ground state energy of the two P donor-pair wavefunc- 
tion is lowered by 7 to 9meV, from that of the isolated P 
atoms. 

Furthermore, in part (c) of these figures, and Fig.[S{b) 
we show the energy level spectrum we calculate for our 
two-electron system, using our extended basis for the sin- 
glet and triplet states. In these plots the difference be- 
tween the first set of excited energy levels cannot be re- 
solved, so we have included an inset which magnifies this 
region. For clarity, we only plot the first eight energy 
levels for both the singlet and triplet two-electron bases. 
An interesting feature of these singlet and triplet energy 
levels is that the four cases of strain give very different 
spectra for the higher energy levels. These plots demon- 
strate that the ground singlet and triplet states are well 
separated from the rest of the higher excited states in 
Hilbert space, for all values of the strain parameter, \- 
This is because the ground singlet and triplet states are 
formed from the symmetric and anti-symmetric combina- 
tions of the single donor ground "Ai" states, which are 
much lower in energy than the next excited single donor 
states, the triplet T 2 and doublet E states. However as \ 
decreases this energy gap becomes smaller as the single 
donor ground state is no longer a symmetric combination 
of the six conduction band states. 

Table [H] shows the difference between the ground 
triplet state and the first excited singlet state. Here the 
first excited singlet state corresponds to two-electron or- 
bitals formed using symmetric combinations of the donor 
electrons at both P donors, and there is negligible contri- 
bution from the doubly occupied orbitals. These results 
are in full accordance with the single donor results re- 
ported earlier in Table |U We can see clearly the trend in 
the plots, that as the strain parameter decreases the en- 
ergy gap becomes smaller. However these energy gaps all 
remain much larger than the exchange coupling between 
the ground singlet and triplet states, and much higher 
than ksT ss O.lmeV, at the cryogenic temperatures re- 
quired for quantum computing. Thus we can consider 
that our targeted Hilbert space, the H-L states, are well 
separated from the rest of the excited Hilbert space. 



TABLE II: Energy gap between the ground triplet state and 
the first excited singlet state, AE (meV). 



and attractive exchange integrals. We see that the attrac- 
tive terms between the exchange charges and the nuclei 
outweigh the repulsive terms, and we obtain a molecular 
binding energy that is deeper than the sum of the en- 
ergy of the two isolated P atoms^ The single donor A\ 
symmetry state wave functions that form the "Heitler- 
London" two-electron ground state, are not spherically 
symmetric, and the electron density for these orbitals are 
more heavily weighted along the co-ordinate axes (see 
Fig. |2| • As a result, even at large interdonor separa- 
tions for 14nm < R < 18nm (along the y direction), the 



X 


R = 14.118nm 


R = 17.376nm 





11.822 


11.808 


-1 


8.283 


8.315 


-5 


4.343 


4.380 


-20 


3.339 


3.378 



Tables IIIII and IIVI list the lowest set of singlet and 
triplet two-electron eigenvalues and eigenvectors for x = 
and R — 7.602 and 14.118nm respectively, and the cor- 
responding single donor basis states contributing. Sim- 
ilarly, table shows the singlet two-electron states for 
X = —20. These tables clearly show the difference in the 
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TABLE III: Singlet and triplet energy levels and correspond- 
ing two-electron eigenstates for R = 7.602nm and \ — 0- 



X = , R = 7.602nm 


Singlet energy levels 


Energy 
(meV) 


Two-electron 
basis states 


One-electron 

bL<Alfc;o lllVUlvfcJU- 


-109.18685 
Doubl 


*43 

y occupied *f /*f 2 


Ax 


-97.77240 | 

Doubly occupied *f /\tf /*f 3 /*25 


Ax/l 2 (X) 

and A x /T 2 {z) 


-97.76809 | *44/*49/*4 6 /*61 

Doubly occupied *f /*f /*f 3 /*f 6 


Ax/T 2 (x) 
and A x /T 2 {z) 


-96.38458 | ¥f 7 /¥&/*w/¥ra 
Doubly occupied *f /*f /*f 6 /*f 7 


A x IE 


-96.32310 |*| 5 /*55/*67/*73/*48 

Doubly occupied *f /*f 4 /*f 2 /*f 


Ax, Ax/E 
and Ax/T 2 (y) 


-95.94663 |*&/*J r /*&/*?T/*?3 
Doubly occupied 

*f/*?/*le/*fr/*li/*f 


and Ax/T 2 (y) 


Triplet energy levels 


Energy 
(meV) 


Two-electron 
basis states 


One-electron 
states involved 


-106.90440 






-97.77305 
Doubl; 


*3 2 /*37 

r occupied /*ffxe 


Ai/T a («) 


-97.76847 | *3 4 /*I 9 
Doubly occupied ^ , 3 r /$ig 


Ai/T a (x) 


-96.50967 | *i 6 /*36/*5s/*6i 
Doubly occupied /^f /^Jg/^Io 


Ax/E 


-96.19400 
Doubl; 


*3 3 /*43 

/ occupied ^2/^17 


Ax/T 2 (y) 


-95.76107 | 
Doubly occupied *2 /*s7*i7 


Ax/T 2 (y) 
and Ai/iJ 



eigenvector basis components with and without strain 
applied. For R — 14.118nm, the ground singlet state is 
the Heitler-London state ^ff 3 , which is composed of the 
single donor ground states at Q\ and Qi- But we see in 
table IVI for \ — ~ 20 that this single donor ground state 
is no longer the six- valley degenerate Ai symmetry state, 
as the strain has broken the degeneracy of the six valleys. 
The lowest energy states are when the effective Bohr ra- 
dius in the direction parallel to the strain is decreased, 
ic. the single donor F{ (15) basis states. 

For x = the first two singlet and triplet excited states 
are nearly degenerate and these two-electron eigenstates 
involve significant contributions from the single donor Ai 
and T2 symmetry states. Here we see in Table IIIII and 
IIVI that for x = these two states involve the T2 states 
in the ±cc and ±z valleys, and are lower in energy than 
the next states involving the T2 states in the ±y val- 
leys. This is because the overall two-electron/two- nuclei 
coupling leads to a more stable configuration when the 
donor electron densities are centered toward each other 
along the inter donor axis (y-axis), ie. the ?2(x) and T2(z) 



TABLE IV: Singlet and triplet energy levels and correspond- 
ing two-electron eigenstates for R = 14.118nm and \ — 0- 



X = , R = 14.118nm 


Singlet energy levels 


Energy 


Two-electron 


One-electron 


(meV) 


basis states 


states involved 


_qq Q8734 


* 43 




-88.14626 
-88.14312 




Ai/T 2 (a:) 
Ax/T 2 {z) 


-88.12780 
-88.12567 




Ax/T 2 (y) 


-88.11610 
-88.11399 


*f*/*S>/*&/*fi 


Ax/T 2 {z) 
and Ai/Ta (x) 


-86.83179 
-86.82667 
-86.82358 
-86.80250 


*47/*48/*67/*73 

*f 7 /*«/*67 
*f7/*«/*67/*ra 


Ax/E 


Triplet energy levels 


Energy 


Two-electron 


One-electron 


(meV) 


basis states 


states involved 


-99.96806 


*31 


Ax 


-88.14459 
-88.14316 


*3 4 /#49 
*32/*37 


Ax/T 2 (x) 
Ax/T 2 (z) 


-88.12656 
-88.12504 




Ax /T 2 (y) 


-88.11693 
-88.11485 


*32/*37 


Ax/T 2 (x) 
Ax/T 2 (z) 


-86.83227 
-86.81847 
-86.79819 
-86.79310 


*3 5 /*36/*5 5 /*61 
*M/*M/*«/*81 
*I 5 /*L/*J 5 /*6X 
*3 r 5/*3 P 6/*i 5 /*I'l 


Ax/E 



TABLE V: Singlet energy levels and corresponding two- 
electron eigenstates for R — 14.118nm and x = —20. 



X = -20 , R 


= 14.118nm 




Energy 
(meV) 


Two-electron 
basis states 


One-electron 
states involved 


-255.87009 


^43 




'F^ilS) + F^ z) (lS)] 


/y/2 


-252.49159 
-252.45132 


*f 4 /*f 9 
*f 4 /*f 9 


«F! W (15) 


-249.11343 






'Fl z) {lS)-Ft z HlS)^ 


/V2 


-225.97320 
-225.82787 






Doubly occupied 
V 1 (z) (15)+F 1 ( - Z) (15)' 


/2 



states. In contrast, Table Ivl shows that for \ — —20 that 
all the higher excited states involve significant contribu- 

( is) 

tions from the F{ (15) basis states. This is because 
with such a large strain applied, the lowest single-donor 
eigenstates are when the effective Bohr radius in the di- 
rection parallel to the strain is decreased, ie. the single 
donor F{ (15) basis states. 
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Tables IIIII and IIVI compare the degeneracy lifting of 
the T and E states for small and large inter donor sep- 
arations. This may be interesting for or relevant to the 
proposed optical Raman experiments in the unstrained 
caseiSS We find that lowest few energy excited states are 
formed when A\ mixes with either T2 or E symmetry 
states. For R = 7.602nm, the two-electron A\jE state 
is even lower in energy than the A\[T%{y) state. As we 
noted earlier this is because the T^(y) state has a smaller 
effective Bohr radius along the inter donor axis. The 
two-electron donor pair is more stable when the electron 
wave functions are more centered toward each other along 
the inter donor axis. For example, at the energies corre- 
sponding to the states containing the E one-electron sym- 
metry states, the two-electron wave-functions are very 
mixed and we are no longer able to assign the energy 
levels to pure basis states. 

For the smaller donor separation with % = 0, and with 
a strain applied, the two-electron wave functions con- 
tain significant contributions from the doubly occupied 
orbitals. This study shows that these doubly occupied 
states are important basis functions to include. We have 
improved upon a previous studySsL where only H-L type 
orbitals were considered. Because we use an extended 
basis which includes both doubly occupied and H-L type 
two-electron orbitals, we find that the two-electron wave 
functions are often mixed states. 

The two-electron energies given in Table IVl and plotted 
in Fig. [7| to are not scaled so that the conduction band 
bottom is at zero-energy. We evaluated relative energy 
shifts using the valley strain parameter, x, for the single- 
donor strain Hamiltonian matrix and neglected any shift 
proportional to identity in it, to be consistent with the 
calculations of Koiller et aliii Therefore, the calculated 
energies do not refer to the zero-energy to be at the bot- 
tom of the conduction band. But the energy eigenvalues 
shown in Table IVl and Fig.[7|to[5]give the correct relative 
splitting among the eigenstates. 

Even just considering the four H-M states in our two- 
electron basis, gives an exchange energy which is very 
close to the extended basis calculation. For the range of 
device parameters we consider, the H-M calculation is the 
most convenient method to use since it is relatively in- 
expensive and very accurate. For this reason we use this 
method exclusively in the next section to obtain accurate 
results expediently and rapidly. 



B. Results using Hund-Mulliken basis 

In figurc HUT a) we show the H-M calculation for the ex- 
change energy for a range of inter donor separations along 
the [010] direction to compare with the previous section. 
This plot demonstrates the oscillations in the exchange 
energy due to the inter valley interference between the 
degenerate conduction band minima. This oscillatory na- 
ture of the exchange coupling has already been reported 
by Wellard et aim and Koiller et a h 11 ! 12 using a H-L cal- 
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FIG. 10: Comparison of the Hund-Mulliken exchange cou- 
pling for different values of strain parameter, \, f° r R m the 
[010] or y direction. We also plot the strain dependence of 
the double occupancy probability in (b). 



culation. We have improved upon and checked the H-L 
approximation, by extending our basis to include the H- 
M states, and found that the H-M basis offers some im- 
provement for the close inter donor separations, and also 
allows us to calculate the ground state double occupancy 
probability of both electrons on the same donor. 

We observe that for the zero strain case the oscilla- 
tions in the exchange are the most conspicuous, as the in- 
ter valley interference is highest, because there are equal 
contributions from all six valleys in the ground singlet 
and triplet states. However, when a strain is applied in 
the z-direction, we find that the ±z valleys are favored 
energetically (see Fig. [3]), which implies a larger effective 
Bohr radius along the inter donor axis (y) , and hence the 
exchange coupling is improved substantially for this par- 
ticular orientation of the P donor atoms along the [010] 
axisiii 

In addition to calculating the variation of the exchange 
coupling with inter donor separation in Fig. 1101 we also 
calculated the probability of the ground state double 
occupation of the H-M singlet states (\&f and ^22) m 
(b)ii£ The exchange coupling with only the H-M states, 
matches the exchange coupling with the full spectrum of 
states very accurately. Thus, based on this we expect 
that the most appreciable contribution of the doubly oc- 
cupied states will be from the H-M doubly occupied sin- 
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glet states in our basis, which is why we only report the 
double occupancy probability of these states. 

This double occupation probability is also an impor- 
tant parameter for quantum gate operations using elec- 
tron spins as qubits. For non-zero exchange coupling the 
spin degrees of freedom are also correlated with the or- 
bital motion, and there is some probability that both 
electrons will be on the same donor These doubly oc- 
cupied states are not part of the targeted computational 
space, and as such, are a potential source of leakage er- 
ror in solid-state quantum computers. These plots show 
that as the strain parameter decreases, (and the effective 
Bohr radii in the x and y-directions increases), this dou- 
ble occupancy probability increases as one may expect. 

Schlicmann et a/«l& showed that gate operations on 
coupled quantum dot pairs which temporarily increase 
the exchange splitting, in order to swap electronic spins, 
inevitably lead to a finite double occupancy probability 
for both dots. However they showed that this double 
occupancy amplitude does not lead to significant errors 
in quantum computing, provided that after the gate ac- 
tion is completed, the double occupancy probability is 
vanishingly small. But if the double occupancy probabil- 
ity occurs to any sizable extent before or as a result of 
the gating action, then any quantum computer based on 
this hardware is likely to fail. Thus for the fabrication of 
these devices we want to minimise the double occupation 
probability for all states at zero voltage. 

The exchange coupling increases consequently with a 
uniaxial strain applied in a direction perpendicular to 
the inter donor axis, which achieves faster gating times. 
However, the double occupancy probability also increases 
correspondingly, which increases the error requirement 
subsequently. In the last section we saw that both these 
features also lead to a greater mixing of the ground state 
with the higher excited states, which causes the energy 
levels to become closer together. This energy splitting 
between the targeted H-L orbitals and the rest of the 
excited two-electron states informs us if during the gating 
action, the coupled donor system is well isolated and the 
higher excited states can be safely neglected. This can 
also give us an estimate for gating times, so that the 
gating operation remains adiabatic. 

In Fig. ^2 we plot the four two-electron energy lev- 
els predicted using our H-M basis, in order to compare 
with our more rigorous evaluation of the higher excited 
states in Fig. and El parts (b) and (c). The insets in 
these plots magnify the splitting between the energy lev- 
els. The bottom inset in both plots corresponds to the 
ground state singlet and triplet energy levels, and com- 
pares favorably with the results shown for this splitting 
using the molecular orbital basis, in Fig. E{b) and E{b) 
for R « 14nm. 

We observe that although the H-M basis provides an 
adequate description of the ground singlet and triplet 
"H-L" states and exchange coupling, it is unable to pre- 
dict the higher energy levels accurately. This is because 
as we reported earlier in Table llVl using our full molec- 
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FIG. 11: Plot of the three singlet and one triplet energy levels 
calculated with the H-M basis, for \ = in (a) and \ = — 20 
in (b). 



ular orbital calculations, the first excited singlet state 
does not include significant contributions from the dou- 
bly occupied H-M singlet basis states. This leads to an 
erroneously large energy splitting between the targeted 
ground "H-L" orbitals and the higher excited states. 

Furthermore we also investigated the variation of the 
exchange coupling for Q\ and Q2 displaced at small dis- 
tances away from a targeted inter donor separation along 
the [010] or y-axis. In these calculations we fixed the 
magnitude of the inter donor separation to be 14nm, 
and varied the two angular variables, 9 and 4> defined 
in Fig.H21 

The results presented earlier in this section show the 
exchange coupling for varying magnitudes of R only 
along the y axis. When the inter donor separation con- 
tains non-zero R Xl R y and R z terms, we see a marked 
difference from the relative smoothness of the exchange 
coupling curves with only a non-zero R y term. 

We show the exchange coupling for \R\ = 14nm and 
X = and —20, for varying 9 and <p in Fig. The inter 
valley interference causing the wild oscillations in J(R) 
is highest when the inter donor separation contains all 
non-zero R x , R y and R z terms (i.e. in (c) where both 
9 and <p are varied). Similar results have already been 
reportedSiii using H-L theory, and we have confirmed 
these results using our more extensive H-M basis. 

For x = —20 we see that the exchange coupling changes 
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FIG. 12: We evaluate the exchange coupling for fixed \R\ = 
14nm and small displacements of Q2, about the targeted inter 
donor separation of R = (0, 14nm, 0), by varying 9 and (j> in 
our calculations. 



dramatically for non-zero R z parts in Fig. 1131 (i.e. when 
cf> varies). This is because for \ = ~ 20 the inter val- 
ley interference terms come only from the ±2 valleys, 
as the single donor ground state orbitals favour the F± z 
valleys iii 

If we study the exchange coupling closely we can iden- 
tify where the peaks and troughs occur in the exchange 
coupling as a function of 9 or <fi. For example, if we ex- 
amine Fig. I13f a) we can identify where the peaks and 
troughs occur for x = 0. In calculating the two-electron 
Hamiltonian matrix and overlap matrix elements, we find 
that several of the integrals involving 'J'q™ ( r — m the 
integrand, have a common factor of e l ( k ^— k ^) R ^ m the 
sum over and k v . 

In Fig. I13t a). we are varying R in the xy-plane, and 
consider R — (R x , R y , 0), where R x = 14 cos 9 and R y — 
14sin#, and 37r/8 < 9 < 5tt/8. We recognise that for the 
range of 9 we consider, R y does not vary as much as R x , 
and R x ranges over both positive and negative values. 
Thus we found that R x is the most significant factor in 
determining the peaks and troughs for this plot. 

In the simplest case when 9 — tt/2 and R x = 0, and 
R = (0, 14, 0), there is a peak in the exchange coupling. 
We find that this is due to the fact that the real part 
of e »(*>M— k u) R j g max i mum value of 1, for 18 out of 
the 36 possible combinations of and h v . Similarly 
we observe peaks in the exchange coupling when R x = 
±ma Jk, where m is any integer, and k = 0.85. If fc M = 
±k x = ±(fc,0,0)27r/a°, we find the real part of e ±ik ' R 
equals 1. Again we find that the real part of e l ( k ^— k ^) R 
is 1 for 18 out of the 36 possible combinations of fc M and 
k v . 

It is difficult to determine the magnitude of the peaks 
in the exchange coupling, because this magnitude also de- 
pends on terms involving e ± lk v R = e ±t(2nkR v )/a W ] J \ C ] 1 
is a complicated function of R y and 9. However in gen- 
eral we observe local maxima at values of 9 for which 
R x = ±ma°/fc, where m is any integer. 

Conversely we find troughs in the exchange coupling 
occurring at values of 9 for which R x = ±ma°/(2fc), 
where m is an odd integer. Here we find that the real 
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FIG. 13: Plot of the H-M exchange coupling for \R\ = 14nm 
and x — and -20. In (a) we calculate J(R) in the a;y-plane, 
where 3tv/S < 9 < 5n/8. In (b) we calculate J(R) in the yz- 
plane, where 37r/8 < </> < 5n/8. In (c) we calculate J(R) in 
the [lll]-plane, where 37r/8 <9,<fi< 5n/8. Here the marked 
points correspond to actual data points evaluated using the 
H-M method, the lines drawn are a guide for the reader only. 



part of e'Cv-M -R j s \ f or i q our , f the 36 possible com- 
binations of fc M and k v , and —1 for 8 of the remaining 
combinations. Thus at these values of 9 we observe lo- 
cal minima in the exchange coupling. Because we are 
only able to evaluate the exchange coupling at finite grid 
points, the peaks and troughs were best matched to the 
data points available, thus enabling us to identify these 
trends. 

For x = ~ 20 in Fig. EH a ) we nn d that the ex- 
change coupling is relatively constant. This is because 
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for x = — 20 the dominant contributions in the ground 
state come from the two ±fc z valleys, and only very small 
contributions from the other four valleys. Thus when 
R = (RxjRyjO), the real part of e ±tk * R is 1 for all the 
combinations of ±fc z , and since the ground state has its 
largest components only in the ±fc z valleys, the exchange 
coupling is maximised and almost constant, for this ori- 
entation. The small fluctuations in the exchange coupling 
are most likely due to the fact there may be small con- 
tributions from the other four conduction band minima, 
which oscillate as a function of 9 as we saw earlier. In (b) 
and (c) of this figure we observe that when R contains 
non-zero R z part, the exchange coupling may oscillate 
even more wildly for x — ~20, as a result of the inter 
valley interference from only the two conduction band 
minima, ±k z . 



V. CONCLUSIONS AND FUTURE 
DIRECTIONS 

It is imperative to know precisely, the single electron 
wave function and two-electron states to determine accu- 
rately the parameter regime necessary for a nuclear spin 
or electron spin quantum computer. In this paper we 
provide detailed electronic structure calculations using a 
molecular orbital method, for a pair of P donor electrons 
in relaxed and uniaxially strained Si. 

We have determined the excitation spectrum of two 
electrons on P donors in relaxed and strained Si, and 
studied its dependence on donor positioning in the Si 
lattice. In particular, we concentrated on the targeted 
ground state singlet and triplet "H-L" orbitals, and ex- 
amined the isolation of these ground states, from the rest 
of the excited Hilbert space. Furthermore we calculated 
the exchange coupling and double occupancy probability 
as a function of strain and donor position. 

Both the exchange splitting and double occupancy 
probability have a similar dependence on inter donor dis- 
tance and lattice positioning of the P atoms. Thus, a 
compromise is needed to maintain a very small double 
occupancy probability in zero field, while realizing a siz- 
able exchange coupling during a gating action. In unison 
with previous theoretical studies of the exchange coupling 
of P donors in relaxed and strained Sif^iiiiS we found 
that the exchange coupling, (and thus double occupancy 
probability), is extremely sensitive to the relative orien- 
tation of the two P donors in the Si lattice. However, the 
energy level spectrum appears not to be affected by the 
relative orientation of the P donors, as these energies are 
on a much larger scale than the exchange coupling. 

The oscillations in the exchange coupling due to in- 
ter valley interference, have serious implications for any 
quantum computer architecture that relies on the ex- 
change interaction to couple qubits. This sensitivity can 
be reduced in the presence of strain, for displacements of 
the donors within the plane perpendicular to the direc- 
tion of the uniaxial strain. We have also identified the 



values of R that lead to the peaks and troughs in the 
exchange coupling. 

It would be useful to investigate the effect of an ap- 
plied voltage on these oscillations, to examine the ex- 
change coupling during and after a gate operation, to 
determine what device parameters are required for fault- 
tolerant quantum computation. It is also important to 
study if the gating action can be performed adiabatically, 
i.e. if during the evolution of the two-electron system 
there remains a finite gap between the ground and ex- 
cited states^ 

Wellard et alm& have extended these calculations to 
include the voltage dependence of the exchange coupling 
within the Heitler-London framework. We hope to de- 
velop their results further by using our extended basis 
to calculate the voltage dependence of not only the ex- 
change coupling, but also the energy spectrum of the 
two-electron system, and double occupancy probability 
of the ground states. The energy spectrum informs us 
if during the gating action, the coupled donor system is 
well isolated, and the higher excited states can be safely 
neglected. In addition, the double occupancy probability 
gives us an estimate of the error rate. However, including 
the electric field is a very computationally intensive task, 
as it requires evaluating the basis functions on our grid, 
over a much greater range of device parameters, which is 
by far the most time consuming part for our calculations. 

The results presented here using effective mass theory 
provide a solid foundation for future device modelling. 
Ongoing work on this project is focusing on extending 
the multi valley calculations to investigate the effect of 
gate voltage on the oscillations in the exchange coupling. 
One would expect that as the gate voltage is turned on, it 
will become more favorable for the donor wave function 
to distort toward the gate. As a result, the donor wave 
function no longer has equal contributions from all six 
valleys, and the oscillations may smooth out as the inter 
valley interference effects decrease. We will need to im- 
plement the full molecular orbital approach to obtain the 
higher excited two-electron states accurately. Thus using 
this method we gain insight not only into the exchange 
coupling and double occupancy probability, but also the 
conditions required to perform adiabatic gate operations. 
To make these calculations more tractable, we envisage 
that the Hund-Mulliken basis should suffice in calculat- 
ing the exchange coupling over a greater range of device 
parameters, as we have shown here that this is the best 
compromise between accuracy and speed. 
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